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ABSTRACT 

The mass of unresolved young star clusters derived from spectro-photometric data 
may well be off by a factor of 2 or more once the migration of massive stars driven 
by mass segregation is accounted for. We quantify this effect for a large set of clus- 
ter parameters, including variations in the stellar IMF, the intrinsic cluster mass, 
and mean mass density. Gas-dynamical models coupled with the Cambridge stellar 
evolution tracks allow us to derive a scheme to recover the real cluster mass given 
measured half-light radius, one-dimensional velocity dispersion and age. We monitor 
the evolution with time of the ratio of real to apparent mass through the parameter 77. 
When we compute rj for rich star clusters, we find non-monotonic evolution in time 
when the IMF stretches beyond a critical cutoff mass of 25.5 M©. We also monitor 
the rise of color gradients between the inner and outer volume of clusters: we find 
trends in time of the stellar IMF power indices overlapping well with those derived 
for the LMC cluster NGC 1818 at an age of 30 Myr. We argue that the core region 
of massive Antennae clusters should have suffered from much segregation despite their 
low ages. We apply these results to a cluster mass function, and find that the peak of 
the mass distribution would appear to observers shifted to lower masses by as much 
as 0.2 dex. The star formation rate (SFR) derived for the cluster population is then 
underestimated by from 20 to 50 per cent. 

Key words: methods; numerical - stars: evolution - stars: kinematics - stars: lumi- 
nosity function, mass function - galaxies: star clusters - galaxies: clusters: individual: 
NGC 1818 



1 INTRODUCTION 

Star clusters are traditionally thought of as old primordial 
structures with ages ranging up to a Hubble time, hence the 
emphasis of theoretical modeling on their long-term evolu- 
tion (e.g. Spitzer 1987; see Meylan and Heggie 1997 for a 
review) . However the wealth of massive young clusters with 
spectroscopic ages of less than 100 Myr observed with the 
HST in interacting galaxies (e.g. the Antenna, M81/82) has 
driven much interest in recent years to the understanding of 
their formation and early evolution. Closer to us, the Large 
Magellanic Cloud hosts a set of young clusters (Elson et al. 
1989; Elson 1991), of which some show signs of primordial 
meiss segregation. The cluster NGC 1818 is one such clus- 
ter where colour gradients are difficult to account for other 
than as a result of their formation history (Hunter et al. 
1997). These issues clearly have bearing on the subsequent 
dynamical evolution and photometric properties of clusters 
and their immediate surroundings (e.g., Lamers et al. 2006). 
To understand very young clusters in quantitative de- 



tail poses a particular challenge to theorists since realistic 
models must account both for the dynamics and the rapidly- 
evolving photometric properties of a young stellar popula- 
tion. In this contribution, we aim to identify evolutionary 
trends that lead potentially to large errors when measur- 
ing the mass of an unresolved cluster with photometry and 
spectroscopy, an effect that bears on all observable cluster 
properties. 

When in dynamical equilibrium, the virial theorem 
gives an exact relation between mass M and mean three- 
dimensional velocity dispersion a: 

\W\ ^ rga^ 

a^ ~ G ' 

where W is the gravitational potential energy, G the gravi- 
tational constant and r^ a radius so defined. All quantities 
entering Eq. (0 must be matched with observables in pro- 
jection. The line-of-sight velocity dispersion ai^^ equals a^ /3 
for an isotropic velocity field, while the gravitational radius 
may be expressed in terms of the projected half-light radius 



M = 



(1) 
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Rhi^ as 



- X - _Rhi, 
2 3 



(2) 



where the numerical factor 5/2 gives a rough conversion to 
a wide-range of clusters fitted with a King mass profile, and 
the factor 4/3 comes from projection on the sky (e.g., Mc- 
Crady et al. 2003; Spitzer 1987 §1.2). Equation Q applies 
when light traces mass throughout the cluster. With this in 
mind we may isolate for M in Eq. Q to obtain 



M = r) 



Rhl 0"lo 

G~ 



(3) 



where the dimensionless parameter rj ~ 10. Several authors 
have used 77 « 10 combined with spectro-photometric data 
to derive M from Eq. (|3J . Such mass estimates can be com- 
pared to masses derived from synthetic stellar populations of 
the same King models fitted to the data to set constraints 
on the stellar IMF. For instance, the stellar population of 
massive Antenna clusters appears to be inconsistent with 
a universal (field) stellar IMF (Mengel et al. 2002; Smith 
& Gallagher 2001). And several clusters in the galaxy M82 
are found to be over-luminous with respect to their mass, 
suggesting a top-heavy stellar IMF in these clusters (Smith 
& Gallagher 2001; McCrady et al. 2003). 

The above studies have taken a fixed value of ry for clus- 
ters of ages up to t = 100 Myr. This simplification, while 
intuitively appealing, was shown recently not to be of uni- 
versal use for massive clusters (Boily et al. 2005). Dense, 
populous clusters will fill the entire range of stellar masses 
drawn from the IMF. This has the effect of dramatically re- 
ducing the mass-segregation time-scale compared with the 
relaxation time-scale and driving heavy (bright) stars to the 
centre of the cluster. The measurements of _Rhi and ctios are 
then biased to values associated with a specific stellar pop- 
ulation, and not the cluster as a whole as assumed in de- 
riving Eq. ^. When the density of the cluster is low (at a 
given number of stars), this bias is reduced and rj remains 
constant over time to a good approximation. Using theoret- 
ical gas-dynamical models, Boily et al. (2005) found a rough 
threshold of mean surface density such that when the initial 
cluster density is (E) ~ 10'' M0.pc~'^ or more, masses de- 
rived assuming constant 77 systematically underestimate the 
real mass by a factor of a few. 

This contribution explores a fuller range of parameters 
and addresses other issues (colour gradients, systematics) 
not covered by Boily et al. (2005). In the next section we 
briefiy recall the dynamical time-scales relevant to the prob- 
lem and show explicitly why a bias should be anticipated 
when deriving the mass of rich, dense clusters. In §3, we 
give details of the numerical approach used to conduct the 
study. In §4, we discuss how the models were analysed and 
quantify numerical and systematic errors. §5 presents the 
results of our survey. In §6 we apply these results to the 
profiling of the stellar mass function and colour gradients 
in clusters. We also explore their implication for a cluster 
mass function, and show that the star formation rate in- 
ferred from cluster populations may be strongly biased to 
lower values. The concluding section introduces a diagram 



^ By convention and when possible, projected quantities are de- 
noted with upper case letters. 



that relates observed cluster properties to their underlying 
potential and draws attention to future developments. 



2 THE DYNAMICS OF MASS SEGREGATION 

Two conditions have to be met for Eq. ^ to be applicable. 
First, all stellar components should be in dynamical equi- 
librium, a sensible assumption whenever the cluster age ex- 
ceeds the virialisation time-scale, i.e., several system crossing 
time tcT, where 



2 rhm/o 



(4) 



with rhm being the spherical half-mass radius. The mechan- 
ics of virialisation leads to equilibrium velocity distribution 
functions independent of stellar masses when all stars have 
the same radial distribution. CoUisional gravitational dy- 
namics, on the other hand, sets a trend towards equiparti- 
tion of kinetic energy between stars of different masses as 
the system evolves. The resulting instability has been stud- 
ied by Spitzer (1969), see also Khalisi et al. (2006) for a re- 
cent work. For two-component systems of individual masses 
mi and 7712, this situation is expressed as 

- mi (Ji = - m2 (72 (5) 

and hence the ratio of squared velocities of the stars equals 
the inverse ratio of their masses: heavier stars have lower 
velocities on the mean, and drop to the centre of the cluster. 
The state of dynamical equilibrium is a good approximation 
to the dynamics only when the migration of the heavy stars 
takes place over long time-scales. 

Secondly, the light should trace the mass so that half- 
light and half-mass radii are identical. When this is not the 
case, the mass M may be derived from either of the two 
relations 



M 



Vo 



Rhm 0"mld 



and M = rj 



-Rhl 0"lo 



(6) 



G ' G 

where a^nid is the mass-weighted velocity dispersion in pro- 
jection along the line of sight, and ai^s its light-weighted 
analog, i.e. the line of sight velocity dispersion most directly 
accessible to observation. 770 ~ 10 is the reference value men- 
tioned already. These relations combine to give 



V = Vo- 



Rhn 



(7) 



and hence rj j^ tjo whenever light and mass follow differ- 
ent runs with R. Since bright stars carry all the light but 
a small fraction of the total mass, we expect 77 > 770 as the 
massive stars migrate to the centre, and both i?hm and (Jmid, 
weighted through the near-static total mass distribution, re- 
main essentially constant. 



2.1 Characteristic time-scales 

The total mass distribution of a star cluster evolves slowly 
over the relaxation time ti of single-population clusters given 
by 



0.138 



r-hn 



1/2 



iV 



(8) 



ln(0.4 7V) 
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which is identical to Meylan & Heggie (1997, §7) once Eq. I^J 
is taken into account where we used 2rhm in the definition 
rather than Tg. A'^ is the number of member stars and the 
ratio Thm/'^g ~ 0.4 for a wide range of model fits to observed 
clusters. 

With iV = 500 000, rhm = 4/3i?hm = 1-3 pc (suggested 
from massive clusters data) and a = \/3(Jios ~ 26 krn.s"^, 
we find a relaxation time 



Owing to mass segregation, we anticipate cti^ > CT2^ 



U ~ 1800 icr ~ 180 Myr 



(9) 



and hence no massive cluster with an age of less than 
100 Myr would be expected to show signs of evolution 
due to two-body relaxation. It is this argument that led 
to the widely used assumption of no evolution of clusters 
in young starburst galaxies. However, Farouki & Salpeter 
(1982) pointed out that the trend toward equipartition is 
accelerated as the mass spectrum {nij} of stars is widened; 
their analysis suggests that the mass segregation will take 
place on a time-scale tms given by (Spitzer 1987) 

3/2 

(10) 



(m) p 



ms TT 



rh, 



where p = (M/2)/(47r r^^^-^-^/3) is the mean density inside the 
three-dimensional half-mass radius (an over-line denotes av- 
eraging over space, and brackets averaging by mass); and 
Timax = max{r7ij}, j = 1..J. Note that tms ~ ir when 
the mass spectrum is narrow, i.e., we recover the single- 
component cluster relaxation time. The mean mass (m) 
^ 0.7 M© for a standard Kroupa IMF (Kroupa 2002). Going 
back to the numerical example given in the above, setting 
Timax = 20 M© in Eq. (IIUI already reduces the mass seg- 
regation time-scale fms to a few Myr, suggesting that mass 
segregation will be effective over the life of massive stars. 

2.2 Example: 2-mass-component systems 

Consider a cluster with two stellar masses, mi < 7712, mean 
surface density E and mass weighted squared velocity dis- 
persion (cr^). Each component i has a velocity dispersion 
ai, a surface density Ei and a surface brightness Ai. The 
equilibrium velocity dispersion assuming Eq. J^ leads to 



, 2\ 2 Si / E2 a2 

El ai-' 



E 
2 Si 



, E2 mi 
E v El m2 



(11) 



On the other hand, using light to weigh the quantities yields 

2 2 Ai / A2 mi\ 

c^iw = (^1 -r 1 + 7 ■ (12) 

A \ Ai m2/ 

Since the stellar IMF is peaked at the low-mass end and 
m,i < m2, the surface density of the second component 
E2 < El ~ E. The same quantities weighted by light 
yield a different result. Since mi < m2, the brightness 
Ai ^ A2 ~ A, so that generally A2 m,i > Ai m,2 for a stan- 
dard IMF (see ii4.1.1ll . This gives the following, approximate, 
relations 



(a ) ^ ai and ai„ 



02 . (13) 

By the same line of arguments, we obtain for the radii 

-Rhm ~ -Rhml and _Rhl ^ i?hm2 (14) 



and i?hml > -Rhm2 

gives 



As a result, 77 computed from Eq. (|7J 



V = Vo 



Rhml 0"1 
^hm2 0-2'^ 



> Vo- 



(15) 



This reasoning will hold true for multi-mass cases in 
general since we may replace mi by (m) and m2 by rrimax: 
light- weighted quantities trace positions and velocities of the 
heavier stars, whereas the global potential and kinematics 
are set by the less-massive stars. 



3 NUMERICAL METHOD 

3.1 Gas models 

Large stellar systems share several thermodynamical prop- 
erties with classical gases (Lynden-Bell & Wood 1968). A 
cluster composed of stars of different masses may be likened 
to a set of concentric spheres of ideal gas satisfying Poisson's 
equation. Larson (1970) pioneered a method based on mo- 
ments of the Boltzmann equation by which energy ('heat') 
flows through the system as it would in a fluid. The temper- 
ature of a stellar system, then, is identifled with the local 
square velocity dispersion so that heat may be transported 
from low dispersion regions to high dispersion regions (ow- 
ing to the negative heat capacity of gravity). Stellar colli- 
sions are treated through a local heat conduction equation 
(Lynden-Bell & Eggleton 1980) which may be calibrated 
to give evolutionary tracks virtually indistinguishable from 
those obtained from N-body calculations (see Spurzem & 
Takahashi 1995). 

Bettwieser & Inagaki (1985) give a good insight of the 
hydrodynamical spirit of the model but note that their 
closure equation requires modification for agreement with 
Fokker-Planck models (Spurzem & Takahashi 1995). A com- 
plete and anisotropic formulation based on moments of the 
Boltzmann equation can be found in Louis & Spurzem 
(1991). 



3.2 Integration code 

The numerical code Spedi^ that we use is based largely on 
the formulation for anisotropic stellar systems due to Louis 
& Spurzem (1991). It was developed further by Spurzem 
& Takahashi (1995). The equations are set on a logarith- 
mic mesh using a scheme which is forward-differencing in 
space and centered in time. Time-integration was performed 
iteratively using a semi-implicit Newton-Raphson-Henyey 
method. The gravitational potential is evaluated from the 
updated (total) density profile directly from Poisson's equa- 
tion. 

Spedi has been adapted by one of us (Deiters 2001) 
to include a model of stellar evolution. We refer to the re- 
sulting code as GasTel. Stars are evolved according to the 
Cambridge stellar evolution tracks, which are available in 
a convenient analytical form (Pols et al. 1998, Hurley et 
al. 2000). By the end of their lives, stars have lost a sig- 
nificant fraction of their mass. This mass lost by stars is 



Further details at http : //www . ar i . uni-heidelberg . de/gaseous -model/ 
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Figure 1. Run of half-mass radius for a three-component Plummer model versus time. The three stellar masses 2/5, 1 and 5/2 were 
drawn from a Salpeter IMF. The symbols are for data points lifted from Fig. 1 of Spitzer & ShuU (1975). The dash is an exponential 
decay oc exp (— t/1.57tr). 



expelled instantaneously from the cluster. However, we may 
still compute rj in the approximation that the total clus- 
ter mass remains constant for short evolution times since in 
reahty the gas will not leave the cluster instantly (see b|().3|| . 
All variables are evaluated on a 200-point grid in Heggie 
& Mathieu's (1986) Nbody units. The constant logarithmic 
width between two grid points is dlnr « 0.095. Spatial res- 
olution in the centre is excellent (152 mesh points to the 
initial half-mass radius at 0.6 numerical units) and the grid 
extends up to 60 numerical units. 



3.3 Calibration, tests 

The only free parameter in the equations of the gaseous 
model is the value of the conductivity (sometimes denoted 
A). It is then adjusted to be consistent with N-body calcu- 
lations and to recover the core collapse time in the case of 
a system of A'^ identical masses (Bettwieser & Inagaki 1985, 
§2.2, see also Spurzem 1992). Several tests have been con- 
ducted to compare gas cluster models with Fokker-Planck 
integrations and direct N-body calculations (Giersz & Heg- 
gie 1994; Giersz & Spurzem 1994 ; Spurzem & Takahashi 
1995). Most of these tests were for two-component models 
of modest mass ratios, while the stellar mass range of in- 
terest here covers nearly two decades. We therefore checked 
explicitly that the numerical setup correctly reproduces the 
dynamics of multi-mass models with a broad mass spectrum. 
Spitzer & ShuU (1975) presented results of mass segregation 
from Fokker-Planck calculations of three-component Plum- 
mer models. The three masses were in the ratio 2/5:1:5/2 
and drawn from a Salpeter IMF. Fig.^graphs the time evo- 
lution of the half-mass radius of each component as obtained 
with GasTel (dotted curves) along with the results read off 
Fig. 1 of Spitzer & ShuU (1975). We find very good agree- 
ment with their data. In particular we find the evolution for 
the most massive stars well recovered from an exponential 
decay of the form i?hm(i) ~ -Rhm(O) exp(— i/1.57trh). 



4 METHOD OF ANALYSIS AND ERROR 

ESTIMATES 

GasTel computes for each dynamical mass group the mass 
distribution and velocity dispersion on a radial grid. The 
mass density p{r) and velocity dispersion a{r) of each group 
are integrated along the line of sight to obtain the pro- 
jected distributions at cylindrical radius R using the relation 



(16) 
(17) 



E(i?) = / p{r)dz 

■J —oo 

p{r)a\r)dz. 

-oo 



The half-mass radius and mean velocity dispersion are com- 
puted for each mass group. System averages are then com- 
puted by summing over all groups using either the density 
or the fight fiux as statistical weight. A stand-alone pro- 
gramme is used to pick the most luminous stellar mass at 
each output time. This approach allows us to combine the 
properties of different stars as desired in the analysis, with- 
out having to re-run the simulation with a different setup. 
For example, it will prove illuminating in the first instance to 
monitor physical quantities attached to the most luminous 
stars alone as function of time, before profiling the system 
including contributions from all the stars (see ^4.31 . 



4.1 Mass sampling 

Our overall goal is to describe accurately the early evolu- 
tion of a cluster. This suggests that we seek out a relation 
between the spectrum of stellar masses and the time-scale 
for dynamical evolution given by Eq. IIIOI . before selecting 
a set of stellar mass groups. 
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4.1.1 stellar IMF 

The field stellar IMF in the solar neighbourhood sets a stan- 
dard reference (Kroupa 2002). This distribution function of 
single stars of mass m is well fitted by a piece-wise power 
law, 



fir. 



if m < 1 M0 

if 1 Mq < m < 10 Mq 

if m > 10 Mq 



(18) 



where a = 1.30, P = 2.35, and 7 = 4.0. The value of a has 
significant uncertainties ±0.7 (Kroupa 2002) whose impli- 
cations will be discussed in H5.4I Stellar demographics are 
computed by integrating f{m) dm up from a lower value 
which we set above the brown dwarf limit at 0.10 M©. If 
the real mass distribution extended to M© with the same 
power law, our cut at 0.1 Mq would allow us to account for 
80 per cent of the actual mass and 99.9 per cent of the actual 
emitted light of the low mass stars (M < 1 Mq). The mean 
stellar mass computed from Eq. 1181 is 



{m) = 



/o°°i M0 /('^) dn^ 



0.7 M, 



0- 



(19) 



However, note that due to discrete and non uniform sam- 
pling, in most of the simulations we have (m) « 0.85 Mq. 



4.1.2 The mass spectrum and the importance of stellar 
evolution 

The mean value given by Eq. 1191 is only weakly depen- 
dent on the upper bound of integration. That upper bound 
should be chosen so as to reflect the richness of stellar pop- 
ulations of massive clusters, yet without overwhelming the 
computational scheme. 

The lifetime of a star is a steep function of its mass. We 
find the following polynomial to give a good fit to stellar 
lifetimes in the mass range [5 Mq, 70 Mq]: 



log fiifc = ci X (log(m)) +C2 X log m + C3 



(20) 



where ci ~ 0.96, C2 ~ —3.7 and C3 ~ 4.2; m is expressed in 
solar masses, and tute in Myr. A star will take full part in the 
time-evolution of the cluster through two-body scattering if 
its lifetime exceeds the mass-segregation time of Eq. (IIOII . 
which we rewrite as (dropping the subscript max) 

tms = — (21) 



where 



K 



3 p 



[m) 



rhn 



3/2 



(22) 



Combining the two time-scales allows us to find a reference 
mass to satisfy the relation log tms < logtufc . Substituting 
tiife from Eq. I|20|l we obtain 



< ci X (log(m))^ + (c2 + 1) X log m -f C3 — log K , 



(23) 



a quadratic inequality for log m. Solving for the roots of this 
quadratic we obtain 



2 ci 2 ci 
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The interpretation of this result is straightforward. All stars 
with initial mass m G [m~,m'^] will not migrate much to 
the centre of the cluster in the course of their lifetime on 
the main sequence. Those with masses above m"^ and below 
m~ , will. Therefore to model accurately the very early stages 
of clusters we should ideally include all stars above m+. 
Recently, Figer (2005) has argued from Arches cluster data 
that all stars have initially a mass < 150 Mq: this would set 
an absolute upper limit on the mass spectrum. However the 
impact of such very massive stars on the dynamics is small 
since they carry a minute fraction of the total mass and 
luminosity of the system. Thus, the most massive stars we 
have included in some of the calculations in this paper had 
m = 70 Mq which already exceeds the mass of Wolf-Rayet 
stars. 

The relation of m to cluster parameters is summed 
up in the constant K: the fraction of all the stars that will 
contribute more effectively to mass segregation is, therefore, 
an implicit function of the cluster we wish to model. The 
minimum of Eq. 12311 occurs for 

logm = -^l±i~i.41, (25) 

2ci 

or m ~ 25.5 Mq . This is the only root to the quadratic when 



K = C3 



(C2 + 1)^ 

4ci 



200.0 = K, 



(26) 



where the numerical value follows from our choice of fitting 
parameters but is otherwise uniquely defined. 

The meaning of Kc becomes clear if we recall the def- 
inition of K and t^, Eq. 12111 and JSJ. Note first that there 
are no real roots to Eq. 1231 when K < Kc. Whenever that 
is the case, all stars drawn from the IMF effectively segre- 
gate while on the main sequence and lose very little mass 
in the process. When K exceeds K^, all stars in the inter- 
val [m~ , m+] must evolve significantly on their way to the 
centre. 

But since K increases with the relaxation time ti , itself a 
rising function of the number of stars N (at a given crossing 
time), we may work out a value for A'' beyond which it is 
unrealistic to neglect stellar evolution. We find after some 
algebra that the condition K > Kc reduces to 

p Mq IMyr ln0.4A - °" <= ^ '' 

where we have substituted the numerical factor rhm/rg = 
0.4. Any multi-mass cluster model (N-body or otherwise) 
that does satisfy this inequality and neglects the stellar evo- 
lution processes is in error. We can isolate for N in Eq. 1271 
by taking characteristic values for the crossing time and 
mean stellar mass to be t^r ~ 0.5 Myr and (m) = 0.7 Mq. 
If the stars are not segregated by mass initially then on av- 
erage p/ p = 1 by definition. With these values inserted in 
Eq. lITTl we find A^ > 6 x 10® = Ac, above the census of an 
average star cluster in the Galaxy ((M) = 300 000 Mq), but 
not atypical of clusters in the Antennae (Mengel et al. 2002, 
table 3). 

Note that for a given crossing time, a multi-mass calcu- 
lation with N < Ac will correctly reproduce the migration 
of stars up to a time of order tma even without accounting 
for stellar evolution. Portegies Zwart & McMillan (2002) and 
Giirkan et al. (2004) used this argument to model runaway 
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collapse of a cluster leading to the formation of an interme- 
diate mass black hole. To perform their simulations, they 
supposed that the relaxation time of their cluster was less 
than 30 Myr so that collapse occurs before the first stars 
explode around 3 Myr. 

4.1.3 The choice of mass bins 

To include very-low mass stars in the computation is costly 
and brings little in terms of the time- variation of the light 
curves, cf. Eq. 1101 . The difficulty resides in having to re- 
solve the light profile of the high-mass stars, and the po- 
tential, dominated by sub-solar mass stars, simultaneously. 
We mitigated this problem partly by selecting a low-mass 
cutoff of 0.1 M0, well below the mean mass of Eq. H19|l . We 
then defined the mass ensemble {rrij}, j — 1,2,.. J, such 
that the lifetimes of two successive high-mass components, 
mj,mj+i, differ by ~ 5 Myr using Eq. I)20|l . We refer to 
all stars with initial mass m > 5 M0 as 'massive stars'. In 
contrast, the sampling of the mass spectrum in the interval 
[0.2, 5] Mq followed a geometric mass-doubling progression: 
mi = 0.2 M0, m2 = 0.4 Mq, ms = 0.8 Mq, etc. We define 
mi/2 = mi/2 = 0.1 Mq. In short, the ensemble {rrij} spans 
the mass range 0.2 Mq to mj non-uniformly and allows a 
much-improved focus on the evolution of the massive stars. 
With mass bins so chosen, the IMF is integrated over 



each interval m 



'j+1/2 - 



-1/2 



to distribute the mass within 



each bin j and normalised so that 



^j+1/2 



/(m)dm = N 



(28) 



The geometric mean of two successive mass groups has been 
used to define the bounds of integration {mj±i/2) for each 
group (see Fig. I^J. 

Simulations were done with 7, 14 and 35 mass groups to 
investigate variations due to a finer sampling of the stellar 
mass function (see Fig. |3J . The trends in mass segregation 
are robust to decreasing or increasing the number of groups, 
but are more noisy in calculations performed with on the 
order of only a few groups. All results in this article are for 
runs with J = 35 mass groups unless stated otherwise. 

4.2 Radius determinations and Monte-Carlo 
checks 

We may distinguish between the half-mass radius derived for 
the continuous density profile of the gas model, and the same 
radius derived for an N-body rendition of that continuum. 

The half-mass radius for each component of the gas 
model is computed by integrating once over the entire plane 
to obtain the total mass of the group; the grid is then re- 
sampled to identify the radius i?hm enclosing half of the 
mass. A linear interpolation at the grid points bracketing 
Rhm gives accuracy to second order in the grid interval. 
Owing to a very fine meshing up to and beyond the half- 
mass radius, errors on this radius are negligible. In practice 
one would like to know what errors are introduced when a 
finite-N model is projected on the grid and the same radius 
evaluated from star counts. This is particularly important 
when the number of stars of a given mass group is low and 
statistical ffuctuations comparatively large. 



To that end, we performed two sets of Monte-Carlo 
(MC) tests. First, we computed the half-mass radius for an 
ensemble of A'^ stars from the surface density of a Plum- 
mer sphere projected on the sky. We call the result Rmc. 
Looking at the dependence in N of the fluctuations of -Rmc 
around -Rhm, we concluded that they were of a Poissonian 
form, mostly due to the random selection of a star in the 
given density function. For example, with A = 50 stars, the 
dispersion around the mean value is of order 15 per cent, for 
A — 1000, it is 3 per cent. Second, we perform 1000 Monte- 
Carlo realizations of the 500 000 stars reference model (see 
Table0 using the density functions at 10 and 40 Myr to see 
which are the dispersions of the fluctuations on the mean 
value of the hali-light radius. In both case, the dispersion 
was of order 3 per cent assuring that in a real cluster of half 
a million stars the hali-light radius is dominated by the 1000 
or so (i.e. 0.2 per cent) most luminous stars. 



4.3 Predominant group approximation 

Following the latter remark, and in order to have a better 
understanding of what is going on, we decided to restrain 
our measures to the most luminous component whose half- 
mass radius and velocity dispersion are assumed to be the 
measured half-light radius Rhi and velocity dispersion crios 
of the whole cluster. 

Figure 0] graphs rj for each individual group (thin lines) 
and the brightest stars (thick solid line) . The individual thin 
curves all increase from their value at f = 0. The increase 
is steeper for the most massive stars, as expected. At t « 
10 Myr, these stars become supernovas and turn to faint 
stellar remnants thereafter. The rj we compute for the system 
drops sharply to the underlying value given by the new most- 
luminous stellar population. And so on for each subsequent 
episode of mass loss through supernovae events. Note that 
at later times the rj of individual stellar groups decreases as 
a result of signiflcant mass loss. This trend is again driven 
by mass-segregation, when the lighter remnants are expelled 
from the central region by the massive stars. Such remnants 
behave like point-sources of gravity with no further stellar 
evolution. The trend where rj decreases is caused mainly by 
an increase of the half-mass radius, which is more signiflcant 
than variations in velocity dispersion. 

The value of 77 computed from the most luminous com- 
ponent gives an upper limit on the value of r; derived from 
the integrated light of all stars. An illustration based on 
bolometric light will be discussed in § 16.31 (Fig. 1121 . Note 
that the difference between the two values depends on the 
wavelength of observation. At the near-infrared wavelengths 
often used to study highly reddened starburst clusters, both 
values of r; are essentially identical as soon as the most mas- 
sive stars have evolved off the main sequence. Indeed, short- 
lived red supergiants or AGB stars overwhelm other sources 
of light; at each time, their distribution is that of the cur- 
rently most massive objects. 



4.4 Number of components 

The large oscillations seen on Fig. 01 for a model with J = 7 
components suggested to us to aim for a significantly larger 
sampling of the mass function to reduce noise to acceptable 
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Figure 2. Generic form of the stellar IMF, Eq. 1181 . The reference IMF (solid line) corresponds to the power index triplet (a, /3, 7) = 
(1.30,2.35,4.00). Changing a changes the slope on the left side of the figure. The integration algorithm is illustrated for a mass rrij-. the 
integration boundaries are chosen to be the geometric means mj_i/2 = y/fnjZrpmJ and mj^x/2 = ^/rrij rrij^i which are mid-logarithmic 
intervals. 
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Figure 3. Mean value rrij of each mass group j for simulations with 7 (x), 14 (□) and 35 (^) components. The sampling used for the 
500 Myr evolution simulation ('Fig. ll2l is also plotted (+). 



levels. Boily et al. (2005) had found for a different reference 
model that J = 14 already gave enough precision to identify 
global trends. Some uncertainties remain with J = 14 mod- 
els in the later stages of evolution and in particular when 
the evolution time exceeds 100 Myr. This suggested to us to 
increase J to the largest possible value. The computational 
time however is quadratic in J and so after some experi- 
ments, we settled for a compromise value of J = 35 mass 
groups binned inhomogeneously as described in §4.1. The 
difference between J = 14 and J = 35 models lies mainly in 
a much smoother transition at the time when the massive 
stars undergo rapid mass-loss while the dynamics for the 
same mass-component is less affected. 



5 PARAMETER SURVEY 

We now survey different parameter values for A^, mmax, 
-RhmO and (m). We tried in each comparison to maintain 
all but one parameter fixed to the reference model values 
given in Tabled From equations llOII . JSJ and Q coupled 
with Eq. Q to eliminate a, we get 



(m) 



N 



RhmO 



mm^ ln(0.4iV) ^(m) N/Rk^o 

Vim) NRZ,o 
m^ax ln(0.4 iV) ' 



(29) 



The survey will highlight dependencies of fms on each quan- 
tities. In all the graphs that follow the solid line indicates 
the reference model of Table^unless stated otherwise. 
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N 



. . . Adjustable quantities 

Alio "T-min '"max 

[pc] [Mq] [Mq] 



M 
[Mq] 



Derived quantities 

O'losO VO S(0) 



[knis 



[Mq pc- 



[Myr] 



[Myr] 



35 500 000 



0.2 



20.0 



1.30 2.35 4.00 418 000 



15 



8.6 0.67 X 105 



180 



Table 1. Parameters and useful data for the reference Plummer model. J is the total number of groups, A'^ the number of stars, M the 
initial total mass. ilhmOi ""losO ^'^'^ VO are respectively the initial half-mass radius, line-of-sight velocity dispersion and t). S(0) is the 
central surface density of the cluster, tr and tms its relaxation and segregation times. Note that for J = 35, we have rriiM =0.1 Mq and 
"^35+1/2 = 20.6. We also have m^in = mi and rjimax = "i35. The indexes a, /3 and 7 define the IMF (see Eq. )18| 'l. The 8 parameters 
on the left are adjustable in the code whereas the 6 on the right are derived from thcin. 
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Figure 4. Global rj estimated from the brightest of 7 mass groups. Dashed and dotted curves show r] for each group of mass rrtj as 
indicated. The solid line is the visible estimate for r] chosen as the brightest component at any given time. The sudden jumps coincide 
with the stellar life time of individual groups. 



5.1 Particle number 



out major implications if we are concerned with clusters of 
ages < 100 Myr or so. 



We first investigate the behaviour of rj when changing the 
total number of stars. A'' while the initial half-mass radius 
^hmo is kept unchanged. The number A'^ of the survey ranged 
from 4 to 15 X 10^ stars. These values of N bracket the clus- 
ters of mass equal to the mean mass of Milky Way clusters 
(some 300 000 Mq , Meylan & Heggie 1997) and the very rich 
Antennae clusters of more than 10^ Mq (e.g. Mengel et al. 
2002). The models all have identical mass groups and upper 
mass limit. With the main sequence lifetime of 20 Mq stars 
Ri 10 Myr we expect from Eq. H29|l a more segregated pro- 
file and larger rj at that time for smaller- TV systems, and a 
similar trend in evolution thereafter. Fig.|21graphs r; for five 
values of N in the range indicated. It is clear from that figure 
that richer clusters, with a longer segregation time, show a 
less-rapidly changing 77. This situation carries over beyond 
t ~ 10 Myr when the first supernova events occur, as the 
second most heavy stars continue to converge to the centre 
on their own segregation time-scale, also oc v'iV. As a result 
the mass profiles are less segregated when the stars move 
off the main sequence in succession for runs with higher val- 
ues of N. Overall differences in the profiling of 7] at times 
f > 10 Myr remains small: for instance the average slope 
dr]/dt is ^ 1.20 for the TV = 1.5 x lO" model, and « 1.31 for 
the smallest-N model shown here. The differences are with- 



5.2 Non monotonic evolution ? 

Large- TV clusters will host a very rich stellar population and 
very massive stars. These stars have very large luminosity 
but are extremely short-lived; their impact on the value of 
rj should therefore be more significant on short time-scales. 
It is interesting, then, to follow the behaviour of -q for indi- 
vidual components when the mass spectrum includes heavy 
stars easily identifiable from spectroscopy which could be 
taken as tracers for the global dynamics. A tracer might 
be the brightest stellar component at any given time and 
we have seen how 77 can be estimated from the brightest 
component alone (Fig. 2J- Would rj increase monotonically 
in time if such a tracer was used instead of a global value 
obtained from integrated light? Call rjj the value of i] com- 
puted for a single component j = 1, 2. ..J, of mass rrij, and 
main sequence lifetime tutcj ■ The mass of the brightest star 
is a monotonically decreasing function of t and if j is the 
brightest mass group at time t then we would say rj — rij. 

From basic stellar evolution models we have fufcj+i < 
tutcj, and so r; = rjj in the time interval iiifc-,+1 < t < 
fiifcj. Thereafter rj — rjj-i, and so on. We noted that 77 
is a growing function of t/tms while the stars are on the 
main sequence. Hence to find out whether 77 will increase 
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Figure 5. Evolution of 77 for different total number of stars A'^. Note the very similar slopes of the curves at times t > 10 Myr. The 
stellar IMF was truncated at 20 M0. 



or decrease as we switch from rjj — > r/j-i, it is sufficient to 
check whether tute/tma is a decreasing or growing function 
of stellar mass. In §4, we fitted the logarithm of tiifc/ims 
with a quadratic function of logTTi, Eq. 1251 . We found a 
minimum for the quadratic (i.e. logtiifo/tms) at m « 25 Mq. 
Therefore, if the current brightest stars have a mass that is 
larger than 25 Mq, we should find 77^-1 (tnfej-i) < r^jitmcj)- 
On the contrary, if the current brightest mass is < 25 Mq , 
then ??j-i(tiifcj_i) > ??j(iiifcj) and hence ri{t) measured from 
tracers would increase with time. 

We graph on Fig. Ela) the curves oi rj = rjj for different 
cutoffs of the mass function, ranging from 12 to 70 Mq. 
It is clear that very large-mass tracers completely bias the 
value of rj to large values, however they can only do so for 
very short times running up to ~ 10 Myr. We note that r]j 
is non-monotonic with time for a cutoff exceeding 25 Mq , as 
expected. The mass of a star cluster derived from massive 
stellar spectroscopic tracers (> 30 Mq) is off by a factor 
that can exceed ~ 2 for the reference setup (Table 0. Note, 
however, that at t « 10 Myr and later, all the curves fall 
back on the same profile, to within small fluctuations. For 
comparison, we plot on Fig.|SJb) the solid curve of Fig.|SJa) 
along with the evolution of T^ium, the analog of rj computed 
from light- weighted integrated quantities rather than just 
from the most luminous component. In practice, ryium is cal- 
culated using Monte-Carlo representations of the simulated 
clusters. Clearly, the bolometric ryium < rj a,t all times, as an- 
ticipated from H4.3I If a red filter were applied, ryium would 
be weighted predominantly by red giant stars, the bright- 
est population, and hence the gap between the two curves 
would close up. The time derivative at i > 10 Myr is almost 
unchanged, so that the trends in time that we will derive in 
a forthcoming section applies to either r/. 



5.3 Mean density via i?hmo 

Observed clusters in M82 or the Antennae are compact. 
They show averaged surface densities that may exceed the 
reference value ~ 7 x 10* Mq pc~^ that we have adopted. 



Boily et al. (2005) already noted that the evolution of 77 is 
significant only for clusters with surface density exceeding 
~ 10* Mq pc~^ (at constant number of stars). In another 
context it is known that the mass density of galactic nuclei 
may well exceed lO'^ Mq pc"''. 

We have therefore explored the evolution of clus- 
ters with different central surface densities by multiplying 
lengths by a factor chosen to cover more than a decade 
in density. The results are plotted on Fig. Q for five values 
of projected initial half-mass radius, from 0.5 to 2 pc. The 
low-density model with -RhmO = 2 pc has a central density 
~ 2 X 10* Mg/pc'^ and we find an increase in i] of at most 
20 per cent after 50 Myr of evolution; by contrast the model 
with -RhmO = 1/2 pc of central density ~ 3 x 10^ Mq/pc'^ 
shows a dramatic increase of r; by a factor « 30/8.6 = 3.5 in 
just 10 Myr of evolution. The rapid increase of rj is driven 
by the much shorter dynamical time of the compact cluster 
which trickles down to a shorter tme in Eq. 1291 : icr oc R'^'^ 
implies a segregation time 4'^" — 8 times shorter for that 
model compared with that of the low-density run. 



5.4 Stellar IMF: (m) 

There is much on-going debate concerning the universality 
of the stellar IMF. The shape of the IMF will fix the mean 
stellar mass which enters the definition of the segregation 
time in Eq. 11291 . We already noted that stars at the high- 
end (> 10 Mq) of the mass spectrum carry much light indi- 
vidually but unless their numbers are greatly enhanced con- 
tribute a small fraction of the total mass. In our exploration 
of the impact of the shape of the IMF on the dynamics, we 
have therefore kept the index 7 = 4.0 as for the reference 
setup (Table0, and focused instead on the effect of varying 
the low-mass power index a. As most of the cluster's mass 
is in low mass stars, a dominates the mean mass value and 
bears directly on rj. The same mass range and number of 
groups were used in all cases discussed below. 

Figure |5| illustrates the three different IMFs used to 
perform the simulations plotted on Fig. |S| To encompass 
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Figure 6. Evolution of rj for models with different upper-mass cutoff rrtmax- (a) We set rj = rjj of the brightest component at time t 
(see H4.c{l . When the stellar mass function is extended beyond A/ = 25 Mq, rjj is not monotonic and becomes rapidly very large at early 
times, (b) Same as (a) for mmax = 70 Mq (solid line) and rjium computed from summing the light from all the stars at all times for 
the same model (dotted line). Note that the first 'bump' due to an extended mass range is much attenuated when computed with total 
bolometric light curves from MC sampling. The effect of very bright and short lived stellar states (e.g. AGB) cause many oscillations 
due to the random sampling in the profile of r;iu„j. 



the standard errors of ±0.7 (Kroupa 2002), we varied our 
parameter a by ±1. The chosen upper value a = 2.3 = /3 
((m) = 0.47 M0) corresponds to a Salpeter profile; it is 
reasonable to assume that any stellar IMF must flatten out 
at the low mass end to avoid a divergence in mass. When we 
do reduce a to 0.3 {{m} = 1.3 M0), the shift in the evolution 
of rj is at no time as dramatic as the one for the Salpeter 
value. Flattening the IMF below the reference a = 1.3 profile 
has not a significant effect on 77. 



6 OBSERVATIONAL IMPLICATIONS AND 
LONG-TERM EVOLUTION 

6.1 The Stellar Mass Function 

The shape of the stellar mass function might be expected 
to vary with radius as a result of mass segregation. The 



central region is rapidly overstaffed with high mass stars 
while the outer parts are depleted of them. This trend can 
be quantifled through the power indices a, /3 and 7 of the 
mass spectrum, by comparing the mass function inside and 
outside a reference radius. Unfortunately, a cluster evolving 
rapidly in time offers no fixed reference radius. To palliate 
this, we computed the stellar mass function in two concentric 
surface elements bounded by the half-light radius from the 
most massive group. While not specially meaningful, this 
choice offers the advantage of a direct link with an observable 
quantity. 

From the star counts in each of the 35 mass bins, the 
mass function is retrieved by summing all the mass within 
m and m + dm, and dividing by dm to obtain a density. We 
then least-square-fitted power laws in the ranges [0.1; 1 M©], 
[1 Mq; 10 Mq] and [10 Mg; 20 Mq] as in Table[T] Since the 
binning is not at constant width, we worried that the mass 
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Figure 7. Evolution of 77 for models with different projected initial half-mass radius iJhmO- The denser models have smaller radius. Note 
that the low-density model with -RhmO = 2.0 pc shows hardly any evolution. 
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Figure 8. Evolution of r] for the reference model with different low-mass power index a for the stellar IMF (see Fig.|2j. The mean mass 
are (m) = 1.3 Mq, 0.85 M© and 0.43 Mq, respectively for a = 0.3, 1.3 and 2.3. 



discretisation would introduce large errors in the values of 
the power indices retrieved. To check this, we trained our 
algorithm on the known IMF from star counts at i = for 
the reference as well as a coarser binning: the power indices 
a, P, 7 of Table Q were recovered to ±0.01, which we take 
as standard deviations. 

Fig. 1^ compares the initial mass function (solid line) 
with the mass function derived inside (dashed line) and 
outside (dotted line) the half-light radius at three different 
times. With the standard model, the changes in the mass 
function are small, therefore the model cluster in this sec- 
tion was initially twice as concentrated {RhmO = 0.75 pc) as 
the reference Plummer model (-RhmO = 1 pc) but otherwise 
the same. The three curves are trivially identical at f = 
with slopes given by the Kroupa IMF. As time increases, 
the low-mass power index a remains virtually unchanged in 
the outer region (varying from 1.3 to ~ 1.33) but shows a 



noticeable decrease in the inner part of the system, down 
by ~ —0.25 after 75 Myr of evolution. The power index (3 
for the mass range [1 M©; 10 M0] shows the strongest vari- 
ations of all, down from its initial value by as much as —1 
inside the half-light radius, and up by -1-0.4 outside this ra- 
dius. The steeper slope in the outer region is a direct conse- 
quence of the outward migration of light stars initially inside 
the half-light radius, while heavy stars flow in the opposite 
sense. 

For very massive stars, the situation is made slightly 
more complicated by the fact that the life-time of these 
stars is comparable to, or less than, the evolution times dis- 
played on Fig. |5| The vertical straight line on each panel 
indicates the mass for which the life-time equals the time 
displayed. All stars to the right of this line are low-mass 
remnants from e.g. supernovae events: these stars therefore 
do not contribute to the light profile of the cluster. Stars 
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Figure 9. The mass function after 25, 50 and 75 Myr of evolution for a Plummer model of initial projected half-mass radius 0.75 pc 
(other parameters as in Table HI . The vertical line indicates the brightest group of stars on each panel. The mass function has been 
retrieved, first by computing the projected half-mass radius of the brightest component, and then by binning stars inside and outside 
that radius at each t (see text for details). The slopes a and fS were obtained by linear regression of the data on the left of the vertical 
line. 
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initially in this mass range are now contributing a small ad- 
dition to the census at the low-mass end of the distribution. 
The high-mass part of the diagram is therefore completely 
depleted, and has not been fitted. 

These trends with radius are similar to those measured 
in young LMC clusters such as NGC 1818 (Hunter et al. 
1997; de Grijs et al. 2002b, Gouliermis et al. 2004). This 
cluster has an age of ~ 30 Myr falling between the times dis- 
played on Fig. Eland a calculated relaxation time ~ 250 Myr 
assuming a half-mass radius of 2.6 pc and mass of 30 000 M© 
(de Grijs et al. 2002a) with (m) = 0.85 M©. This is longer 
than in our simulations. Note however that the mean density 
of the model with -Rhm ~ 0.75 pc matches well the density 
inside the half-mass radius of the cluster NGC 1818 (Elson 
et al. 1987). Inspection of Fig. 9 of Gouliermis et al. (2004) 
shows that the power indices in the inner part (e.g. < 0'.3) 
derived from their data are similar to those of our simula- 
tion. This raises the possibility that dynamical mass segre- 
gation may yet play a key role at the heart of that cluster 
while primordial segregation is needed to explain the exter- 
nal parts, a conclusion already reached by de Grijs et al. 
(2002b). 

6.2 Colours 

It is interesting to investigate to a fuller extent observable 
consequences of mass segregation. To that end, we extracted 
colours from our model clusters by coupling the Cambridge 
evolution tracks to the spectral library of Lejeune et al. 
(1997, 1998). Nebular gas emission lines, thought present in 
embedded young clusters (Anders & Fritze-v. Alvensleben 
2003), are left out of the current analysis. 

Sampling the mass spectrum requires some care in or- 
der to minimize errors in colour magnitudes as discussed by 
Chariot & Bruzual (1991). First, we tabulated the various 
evolutionary epochs for a large set of masses split in equal 
logarithmic intervals from 2 Mq to 100 M© . The luminosity 
function was constructed by carefully integrating the light 
flux from all the stellar masses in a given evolution phase, 
paying great attention to resolve such brief but very bright 
phases as the upper asymptotic giant branch (AGB). 

Colours were computed in different wavebands (B, V, 
I from Bessell 1990, and K from Bessell & Brett 1988) and 
compared with those of other authors, who used the evo- 
lutionary tracks of the Geneva group, the Pad ova group or 
variations thereof (Girardi & Bertelli 1998, Bruzual & Char- 
lot 2003, Mouhcine & Langon 2003 and references therein; 
mostly based on tracks by Bressan et al. 1993 or Schaller et 
al. 1992). We found that the Cambridge tracks produce sig- 
nificantly redder colours than others, which have in general 
been more specifically tuned to reproduce the observed inte- 
grated colours of star clusters. The origin of these differences 
lies in the time stars spend in the late, red phases of stellar 
evolution. In calculations using the Cambridge tracks, the 
predominance of red supergiants and upper AGB stars at 
near-IR wavelengths is probably exaggerated. To illustrate 
the effect of the luminous red stars, we ran simulations with 
and without stars on the thermally pulsing AGB (TP-AGB). 
Fig. llOl shows the evolution of the integrated V — K colour of 
the model cluster. The two curves bracket values commonly 
found in the literature for simple stellar populations. 

Colour gradients were quantified by measuring the dif- 



ference between the colours measured inside and outside the 
projected half-light radius: 



Av-K = {V~K)^R,,-{V-K)^R, 



(30) 



Ay-K is positive when the inner half of the cluster is redder 
than the outer half. Fig. 1111 shows the evolution of Av-k 
for the reference model (Table 0. At early times, mass seg- 
regation tends to make the centre bluer, as the massive stars 
that fall towards the centre are still on the main sequence. 
Once these stars evolve off the main sequence, the centre of 
the cluster very rapidly becomes the reddest part. Further 
evolution of the colour gradient presents fiuctuations that 
refiect the lifetimes of stars of progressively smaller initial 
masses and the mass-dependent time they spend as luminous 
red objects. However, Av-k remains above 0.05 magnitude 
throughout, even when the brightest stars such as those of 
the TP-AGB are artificially switched off. 

6.3 Long-term evolution 

It is natural to ask whether the continued increase of r; 
observed for the reference model over a time-scale of ~ 
100 Myr carries over to longer time-scales. Recall that the 
relaxation time of that model ~ 200 Myr from Eq. (|HJ. We 
would anticipate some increase in tj for any cluster with an 
age > ir, however longer relaxation times also imply longer 
mass segregation times and significant contributions from 
low-mass stars to rj. Since the stellar winds of lower-mass 
stars moving off the main sequence are less energetic, it is 
not clear whether or not the residual gas will be evacuated 
on a very short time-scale and how this will impact on the 
dynamics. We take the view that much of the gas can either 
remain within the cluster for a period exceeding tr, or on 
the contrary be evacuated rapidly, thus hopefully bracket- 
ing all realistic cases. Below we assess whether either limit 
(or both) will result in a drop in tj over long time-scales. 

We ran the reference calculation (Table for up to 
500 Myr. That stretch of time would correspond to several 
revolutions through the galactic potential of a galaxy such as 
the Milky Way, and overrun the star-burst phase of a typical 
galaxy merger where young clusters are observed to form. 
Since heavy stars play only a minor role beyond ~ 30 Myr 
of evolution (cf. Fig.|SJ, we split the 35 mass bins so that the 
stellar MS lifetimes now differ by ~ 20 Myr from one bin to 
the next (instead of 5 adopted earlier, see Fig. EJ. In this 
way the part of the mass spectrum below 5 Mq is far bet- 
ter sampled than previously. Nevertheless the discreteness 
of the mass function will cause large fluctuations, particu- 
larly noticeable at times t > 200 Myr (cf. H4.3II and tend 
to underestimate the increase of r; as we have seen when 
comparing simulations with fewer components. 

On Figure 1121 we graph rj evaluated in three different 
ways for comparison. The dotted curve (denoted '77 with 
ML') assumes instantaneous evacuation of the mass released 
through stellar evolution. The dashed curve denoted rji^m 
also assumes instantaneous mass evacuation but is the ana- 
log of 77 computed (as in Fig|SJ from light-weighted quanti- 
ties, using a Monte-Carlo representation of the cluster. To 
verify the impact of stellar mass loss, we also plot the value 
of 77 obtained with a constant total cluster mass. The re- 
sult is shown on Figure [re labelled as '77 without ML' (solid 
line). For that curve we find an increase of r; of a factor 
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Figure 10. V — K plotted against time for the reference model cluster (see Table0. The solid line is for integration with all possible 
stellar phases including the TPAGB. The dotted line is the same but without the TPAGB phase. 
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Figure 11. Difference of colour Av_x between internal and external part of a cluster relative to its half-light radius at each given time. 
The inner region becomes bluer over the first 10 Myr of evolution, when the most massive stars in the sample (mmax = 20 Mq) become 
red giants. Thereafter the inner part remains systematically redder than the outer region by as much as 0.2 dex over the first 50 Myr. 
Note however that the color index Ay_x fiuctuates wildly. The trend of increasing Ay_x at 4 > 40 Myr is continued in time exceeding 
60 Myr (off scale). To account for possible bias from the TPAGB phase, when stars are very bright in the red, we recomputed the colors 
by removing all stars in that phase of evolution (dotted line on the figure). The behavior stays unchanged. 



ranging from 5 to 6 by the end of the run. This curve re- 
mains significantly and systematically higher than the other 
two throughout the run. These curves bound, therefore, all 
possible values because in practice some of the gas will leak 
out and r] oc M would decrease as a result. A full simulation 
including hydrodynamical effects would trace rj somewhere 
between the curve shown here for 77ium and the sohd line of 
r] without ML. 

It is highly likely that tidal fields will affect the mor- 
phology of a cluster as it orbits the host galaxy. This will 
certainly change the profiling of r] in time, for example by 
stripping some of the cluster's mass. But unless the cluster 
evolves in a very strong tidal field, the dynamics inside the 
half mass radius should prove relatively robust. If the to- 



tal cluster mass decreases, stripping mass outside rhi, then 
r] would decrease from the values obtained here. A full in- 
spection of this issue would require three-dimensional model 
clusters which lay beyond the scope of the current study. 



7 APPLICATION TO CLUSTER MASS 
FUNCTIONS 

T.l T) — tro relation 

The results of >I5 . 1 1 and 15.31 suggest a common thread link- 
ing models of high mean density and those of smallish total 
mass. Since these two quantities combine to give the sys- 
tem relaxation time in Eq. Q, we may hope to relate r; 



© 2006 RAS, MNRAS OOO.IHOni 



Mass of Dense Star Clusters in Starburst Galaxies 15 



1U 


1 1 1 1 

rj with ML /\ / 


35 


^lum with ML /\ /\ f \r 


30 


7] without ML /\f^f\\ N N- 


25 


A/^/^ .■■■■■. .V 




/\/^ ■■■,.• ,-."~ ■-.'' 


20 


/-^^^/.'■".:-<.- ■■>'---'''"""''""'' "' 




rJ,'''~~''''' ""-'' 


15 


v^'' 


10 


^'" 


5 
n 


1 1 1 1 



100 200 300 

Time [Myr] 



400 



500 



Figure 12. Parameter r) vs time for an 35-component simulation (see text for details of the numerical setup). The top-most solid curve 
was obtained by computing r\ at constant total cluster mass, so disregarding mass loss (ML) due to stellar winds. The dotted curve 
assumes on the contrary that this ML leaves the cluster instantaneously. A more realistic situation where the stellar winds escape over 
a finite time interval would pitch r\ somewhere between the two curves shown here. Note that r] was once more calculated by taking the 
most luminous component as tracer. If we draw Monte Carlo realisations using the luminosity from all the stars to compute i] = J7iuni 
under the same assumption of instantaneous ML, we instead obtain the dashed curve which shadows closely the dotted line. Either way 
of computing rj yields continued increase up to 500 Myr. 



for models of different relaxation times but comparable ages 
for a given IMF through a simple scaling formula involving 
the initial relaxation time fro- Thus we take the view that 
evolution will be driven almost exclusively by two-body re- 
laxation and not e.g. variations in the stellar mass function. 
Such effects would play a significant role, as shown on Fig.|Sl 
but only for clusters not older than a few million years. 

We sought out such a scaling relation by performing a 
series of runs for models with different values of t^o, span- 
ning a wide range of values in -RhmO and A'' but in other 
respects identical to the reference model (Table0. We com- 
puted rj for these models at t = 10 Myr, approximately when 
the first stars become supernova, and a time f = 40 Myr 
which gives an intermediate age between e.g. age estimates 
of Antenna; clusters (Mengel et al. 2002) and that of M82- 
F (Smith & Gallagher 2001, McCrady et al. 2005). Fig. IBl 
plots the variations in 77 at these two times relative to its 
initial value, Arit/rio, as function of tro- Both sets of points 
are well fitted by single power laws. 






At X t,7* 



(31) 



and we list the values of At and power index at in Table El 
for two different times. The error bars shown on the figure 
result from oscillations of rj due to the mass sampling. The 
power-index at shows only a mild dependence on the age 
of the cluster. This is somewhat surprising if we note that 
the same power-law functional fit applied to the half-light 
radius A_Rhi/_Rhio and square velocity dispersion Aai^/ai^Q 
give equally good results but now the parameters vary much 
between the two chosen times (cf. Table |2|l. 

A general expression for rj valid at ages between and 
about 50 Myr can be derived from the above. Guided by the 
Eispect of Figs. 121 17| or |H| we distinguish an early regime of 
rapid evolution, up to ~ 10 Myr, and a subsequent regime 



of slower changes. A good fit to model values is obtained 
by running a straight line from the value of 77 at t = to 
its value at i = 10 Myr for the early regime, and another 
straight line through the values at 10 and 40 Myr for the 
later evolution, using Table [51 It can be summarised, with t 
in Myr, 



at 



1.36 



At = 



50 i 

lOi + 400 



if t < 10 Myr 
if t > 10 Myr. 



(32) 



Fig. 1141 compares the analytical values with the results of 
a simulation with iro — 205 Myr. Differences do not exceed 
10 per cent even when extrapolating to ages of ~ 100 Myr. 
This level of error is found for all models with an initial 
relaxation time above 100 Myr. When the relaxation time is 
shorter, we find that the interpolation scheme systematically 
underestimates rj. 



7.2 Application to a model CMF 

Using the linear interpolation scheme of the preceding para- 
graph we may compute the current mass conversion factor ri 
of a cluster from Eq. (I3H given iro and its age t < 100 Myr. 
The real mass of the cluster so retrieved may then be com- 
pared to the mass estimate we would have computed had we 
kept the initial value 770 ~ 8.6 constant throughout; the ra- 
tio of real to estimated mass equals rj{tro,t)/rio. To compute 
rjitio, i) first requires a cluster mass and half-mass radius, in 
order to evaluate tro from Eq. JSJ. To do this for an ensemble 
of clusters, we set up Gaussian distributions for M and -Rhm 
as well as the cluster age, t. These distributions were each 
sampled independently with 10^ realisations using a stan- 
dard Ulam-von Neumann (Monte Carlo) rejection method. 



© 2006 RAS, MNRAS 000,0021 



16 J.- J. Fleck, C. M. Boily, A. Langon and S. Deiters 

0.8 



0.7 
0.6 
0.5 
A77 0.4 
^ 0.3 
0.2 
0.1 



1 


'Ar/io 


A7]4o ' 




l-H 


1 — — 1 


1 — h 


— 1 




^0 


Vo 




, H-l' 


500 t^o^-^*^ — - 


800t7o^-36 --- 


- 


5 






- 


$ 








\ %. 








\ •* . 






- 


■^5. 
$- 




■•^■•■■■.£. 


- 


1 


1 1 


-■^ $- 


■ - J 


1 


— s 



loo 



200 



300 



400 500 
iro [Myr] 



600 



700 



800 



Figure 13. The parameter 77 as a function of initial relaxation time ty-o from different simulations where both A'^ and -RhmO have been 
varied. The relative increase Arj/rjo is a power law of irO. The error bars are deviations from the mean value of r] at each time. 

t = 10 Myr t = 40 Myr 



At 500 -90 -57 800 -76 -36 
at 1.36 1.15 1.22 1.36 1.05 1.08 



Table 2. Power law fits to 77 as function of the initial relaxation time. A least-square fit of the functional form defined by Eq. 1311 was 
performed at two cluster ages. The standard deviations on At and at are respectively 5 per cent and 1 per cent. 



In the following, we quote the dispersion a of these distri- 
butions as uncertainties on the mean value (i.e., meanzba). 

The results are shown on Fig. 1151 with a log-Gaussian 
CMF of mean mass 5 x 10"" Mq for two realisations: (a) an 
ensemble of compact clusters of mean radius 1 ± 0.2 pc and 
age 30 ± 5 Myr, in line with values adopted for our refer- 
ence model; and (b) an ensemble of mean radius 2 ± 0.4 pc 
and mean age 60± 10 Myr inspired from M82-F cluster data 
(McCrady et al. 2005). The resulting relaxation time distri- 
butions are respectively 190 ± 70 Myr and 550 ± 130 Myr. 

Because the mass derived assuming a constant ri = rjo 
is always lower than the actual cluster mass, the distribu- 
tion shifts to lower masses as compared to the true CMF. 
Table |3 list the displacements of the peak of the CMF for 
four Gaussian distributions of different mean initial relax- 
ation times tro, each of the same 50 Myr dispersion. It is 
clear that the shorter relaxation times lead to a larger shift 
and we find a maximum shift of 0.2 dex for the distribution 
of average frO = 150 Myr. Very similar conclusions apply 
for the ensemble of longer-relaxation time when the average 
age is also longer (case [b] above, displayed on the right-hand 
panel on Fig. llSH . 

In this spirit the very massive and young Antennae clus- 
ters are particularly interesting. Table 3 of Mengel et al. 
(2002) lists parameters for five young clusters of ages rang- 
ing from 6 to 10 Myr and masses from 600 000 to 5 x lO'' Mq. 
The projected half-mass radius of these clusters is signifi- 
cantly larger than 1 pc, the value we have adopted for our 
reference model. Those large radii may mislead one to expect 
much larger relaxation times and, consequently, little or no 
evolution of the 77 over time. However, we note that Mengel 



et al. (2002) fitted King models with a concentration pa- 
rameter ^/(T Ri 6 to the fight profiles of their clusters. Such 
King models are significantly more concentrated than the 
Plummer model that we used in the calculations performed 
here. In fact, it turns out that a King model with '^ /a = 6, 
a half-mass radius of 2 to 3 pc (cf Table 3 of Mengel et al. 
2002) and a mass of some 500 000 M0 , has a mean density 
within its core radius that essentially equals the mean den- 
sity within the half-mass radius of our reference Plummer 
model. Consequently, the rapid evolution of rj found in this 
article should be applicable to the dynamics in the core of 
some Antennae clusters, implying that the core regions of 
very rich clusters are still affected by strong segregation de- 
spite their very low age, a conclusion also reached by de Grijs 
et al. (2005) using observational arguments. As a result, the 
core will appear more compact, while the half-mass radius is 
left largely unchanged, a situation that leads to King model 
fits with higher values of "if /a than is appropriate (cf. also 
Boily et al. 2005 for examples of this effect). 



8 DISCUSSION 

This paper investigated possible biases when estimating the 
dynamical mass of young and dense stellar systems from 
spectro-photometric data. Using the virial theorem, one may 
convert observed hali-hght radius and flux- weighted mean 
velocity dispersion to mass through the dimensionless pa- 
rameter 77 defined in Eq. Q • This factor will vary with time 
due to mass segregation whenever the two-body relaxation 
time in Eq. Q is short: the heavy bright stars segregate 
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Figure 14. Example of the reconstruction of tiie evolution of ri in time. The symbols {'&) represent the data from the model while the 
various curves are drawn from the parameter fits of Table lousing 3 values of t^o in Eq. I31i . 

(Uo) ± 50 Myr 150 250 350 450 



shift of (log A/) 0.2 0.15 0.1 0.05 



Table 3. Shift of a lognormal distribution in mass centered on 10® Mq for different mean relaxation times. The cluster population has 
a Gaussian age distribution of moan 30 it 20 Myr. 



to the centre rapidly. A parameter space exploration led us 
to conclude that for clusters where frO ^ 200 Myr rj may 
increase by a factor of 2 compared with its initial value, 
whereas if tro > 500 Myr then very little evolution of 77 will 
take place within the first 100 Myr. Meanwhile, the mass 
distribution and potential, dominated by fainter stars, re- 
mains largely unchanged, so that hght does not follow mass 
anymore. 

We can synthesise the main features of this bias in 77 
in a diagram of cluster age versus relaxation time derived 
from observations. Substituting the projected half-mass ra- 
dius Rhm and the total mass Af for r^ and A*' in Eq. ||HJ, we 
get 



irO 



■ -fthm C"mld 



G (m) \/3 1nA 



(33) 



where (m) is given by the IMF and In A is the Coulomb 
logarithm (for which we adopted ln[0.4A'^] previously). The 
same formula applied to observed quantities would give an 
estimated relaxation time 



1/2 



CfjObs 



2x 0.138 ?7o 
G (m) \/3 1nA 



-Rhl 0"los. 



(34) 



With the ratio 77/770 given by Eq. Q we obtain a ratio of 
'true' to measured relaxation times, 



I-r 



t^T,ohs 



0"lo 



0"mld 



(35) 



The ratio of squared velocity dispersion was found not to 
decrease by more than 10 per cent throughout the simula- 
tion time; we may therefore set crios/o-mid = 1 in Eq. (I35II . 



Equation 11351 together with Eq. I31II can be solved with in- 
put cluster age and measured relaxation time to obtain a 
unique pair (trOi^?)- It is then straightforward to draw lines 
of constant 77 in a graph of age vs ir,obs axes and recover the 
true relaxation time of the mass profile (since we must have 
tio = tr.obs at time t = by construction and the potential 
does not change). 

Fig. 1161 graphs the contours lines of constant rj/rjo in 
the plane ir,obs-age. All clusters start off on the age = 
axis which coincides with the contour 77/770 = 1. As the clus- 
ter becomes older and mass segregation sets in, the time 
evolution marks a path that is seen to drift to shorter (mea- 
sured) relaxation times and larger 77. Each level is indicated 
on the graph. The path for our reference model is shown 
along with a second model of initial half-mass relaxation 
time of 400 Myr. Note that even for this model tr.obs drops 
to « 200 Myr after 100 Myr of evolution: at that time its 
mass is underestimated by ~ 30 per cent. Fig. llGr b) zooms 
in on the interval [0 : 100] Myr of the fr,obs-axis. 

As seen on Fig. 1161 the most dramatic evolution in 77 
occurs in the first 10 Myr. Furthermore, if only the most 
luminous stars were used as tracers, still higher factors 
77/770 would be expected (cf. Fig. |HJ . For the young massive 
Antennae clusters for which the measured relaxation time 
tr.obs > 500 Myr, mass segregation is negligible when ap- 
plied to these clusters as a whole. We noted, however, that 
our reference model provides a good fit to the core region 
of some of these clusters in terms of mass and density. Thus 
mass segregation may yet prove an efficient agent for evolu- 
tion in the central part of massive young clusters, a process 
that would contribute to make the core look more compact 
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distribution of initial relaxation time are (a) 190 it 70 Myr and (b) 550 it 130 Myr. 



than it really is and so inflate the concentration parameter 
c = log(r-h/r'c) of King model fits to these clusters. Examples 
of this phenomenon are given in Boily et al. (2005). 

Light curves have been used to estimate colour indices 
averaged over two surface elements bounded by the pro- 
jected half-light radius. The difference between these colour 
indices taken as a function of time shows that the inner part 
becomes bluer by 0.05 dex so long as massive stars have not 
reached the red giant state. At that time and for all times 
thereafter, the inner region becomes redder (by more than 
0.05 dex). Colours should be interpreted with caution. The 
evolution ofV — Kis highly sensitive to the properties of the 
red stages of stellar evolution, and especially to the red su- 
pergiant and asymptotic giant branches. Such fine details of 
stellar evolution, often model-dependent and difficult to pin 
down with precision, have less bearing on the colour gradi- 
ents because these giants stars dominate the light in the red 
wherever they are, and hence Av-k quantifies their concen- 
tration in space. The stellar evolution tracks used here tend 
to exaggerate the role of bright red stars. Discarding the 



most luminous red supergiants or TP-AGB stars altogether 
reduces the gradients by a few hundredths of a magnitude 
at most. 

The evolution of the stellar mass function within the 
model clusters was also quantified through variations in the 
power- indices defined for the IMF Eq. I18II . We noted that 
the variations of the power index /3 during evolution are a 
good match to those observed in the LMC cluster NGC 1818 
(Gouliermis et al. 2004). However the initial relaxation time 
of the more concentrated model used in this section of ~ 
115 Myr is significantly smaller than the one derived for this 
cluster (we compute ~ 250 Myr from Eq. (|5J; see also Elson 
et al. 1987) which implies that dynamical mass segregation 
alone does not account for gradients in the stellar population 
for the cluster as a whole. Thus we would argue that a fair 
degree of primordial stellar segregation must be relied on to 
explain that cluster's photometry. Despite this caveat, it is 
well worth repeating that the central relaxation time of NGC 
1818 of some 120 Myr is of order of the half-mass relaxation 
time of the model used, and therefore dynamical segregation 
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Figure 16. Contours of constant rj{t)/rio (dotted curves) in the plane of cluster age vs relaxation time, tr,obs, derived from cluster 
observables. The contours were constructed using the bi-linear interpolation scheme of H7.ll The contour r]{0)/rjo = 1 coincides with the 
horizontal axis (cluster age = 0). (a) The solid lines trace the evolution in time of rj/rjo for two model clusters with initially t^ ot,s = 180 
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now one of the models has t^ ^^^^ = 100 Myr initially. The evolution track for that case crosses contours of yet higher values of r](t)/rio 
at fixed age. Note the change of scale on the absciss. 



surely has played a role in the evolution of the stellar mass 
function near the centre (see also de Grijs et al. 2002b). 
Our view is that a set of models tailored to that particular 
cluster will be required to disentangle fully primordial from 
evolutionary effects. 

With r]{tro,i) derived from Eq. I31II . it was possible to 
construct a Log-Gaussian cluster mass function and carry 
out a statistical survey of the impact of mass segregation on 
the shape of the observed CMF obtained from assuming no 
evolution of 77, to the CMF derived from taking into account 
the time-evolution of ?7 (Fig. 1161 . The actual total integrated 
mass of the CMF is 50 per cent larger than the apparent 
one for the case describe on Fig lldf a) and 20 per cent for 
Fig llSf b). This has direct bearing on the global star forma- 
tion rate (SFR) derived for galaxy mergers and starburst 
galaxies in general. We noted that the two CMF's so con- 
structed differ mostly at the low-mass range (< 500 000 M0) 
where the relaxation time is significantly shorter. To quan- 



tify this effect, we found a shift in the peak of the real CMF 
towards high masses compared to the 'observed' CMF. This 
shift is on the order of 0.2 dex for a distribution of relax- 
ation times centered around 200 Myr and is lower when this 
mean relaxation time is larger (Table El . Furthermore we 
also found a slight widening of the observed CMF, by a 
logarithmic factor ~ 0.05 (see Fig. I15II . This would have 
some influence on evolutionary predictions of cluster mass 
functions. Vesperini (1998) has shown that the Milky Way 
CMF may well shift to smaller masses by ~ 0.1 dex over a 
Hubble time due to tidal destruction and other effects. The 
trend we found here goes in the opposite direction, however 
it is only operative for clusters with short relaxation times. 
As clusters are possibly more massive than estimated from 
observations but also less concentrated (lower King fitting 
parameter) it is not clear how tides and other disruptive ef- 
fects will shape up the CMF, especially if the host galaxy 
itself is out of equilibrium. A set of fully three-dimensional 
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N-body simulations could enhance our knowledge of the in- 
fluence of mass segregation on longer time-scales and with 
strong tidal fields found for example in merging galaxies. 

Our models of isolated clusters suffer a few important 
limitations. We have mentioned the role that tidal fields will 
play in removing weakly bound stars. Another aspect of the 
problem is the possibly low star formation efficiency when 
the cluster forms. We mentioned how gas from stellar winds 
might impact on the dynamics ( Hfci.^ll . Residual gas from 
the formation epoch will also drive much evolution in the 
early stages by bringing the cluster out of virial equilibrium 
(Elson et al. 1989, Kroupa & Boily 2002, see also Bastian & 
Goodwin 2006). Yet another aspect is our tacit assumption 
that stars are all born at the same time and all evolve in 
unisson. Stars in massive clusters may well have ages that 
vary by a few million years. This will have some bearing 
on the rise of rj in the early stages because not all the stars 
become remnant at the same time, leading to enhanced mass 
segregation and further increase in rj. For instance, the knee 
seen at t ~ 10 Myr may well increase to yet higher values 
before shifting over to the slower rate of increase that we 
have advertised (Fig.|SJ. 

A more severe limitation however, one that will impact 
on ri at all times, is the fraction of primordial binaries. Tight 
binaries will survive for eons and in particular a very large 
fraction of them will survive for the short times that are 
of interest here. The presence of binaries and multiple stars 
naturally enhances the observed velocity dispersion which 
biases the mass estimate to larger values through the virial 
theorem. However, binaries also instantly broaden the width 
of the effective stellar IMF, since, roughly speaking, they will 
dynamically act as single stars of mass equal to the sum of 
their member stars. If the fraction of primordial binaries is 
low, the mean stellar mass will remain unchanged but the 
maximum mass will effectively double. The net effect, then, 
is similar to halving the mass segregation time-scale ims by 
reducing the mean stellar mass. This can be accomplished 
by increasing the power index a of the IMF in Eq. 11811 . We 
have found after ~ 10 Myr of evolution for the extreme case 
where a = 2.3 (Salpeter value) that rj has nearly increased 
by 30 to 50 per cent in comparison with the result for the 
standard case a = 1.3. These considerations clearly point 
to yet more rapid evolution and the need for more realistic 
models than is affordable here to pin down more precisely 
the dynamics of young massive clusters. 
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